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Abstract 

We study the spectral statistics of the Dirac operator on a rose-shaped graph — 
a graph with a single vertex and all bonds connected at both ends to the vertex. 
We formulate a secular equation that generically determines the eigenvalues of the 
Dirac rose graph, which is seen to generalise the secular equation for a star graph 
with Neumann boundary conditions. We derive approximations to the spectral pair 
correlation function at large and small values of spectral spacings, in the limit as 
the number of bonds approaches infinity, and compare these predictions with results 
of numerical calculations. Our results represent the first example of intermediate 
statistics from the symplectic symmetry class. 

1 Introduction 

For classically chaotic quantum systems with an anti-unitary symmetry (for example, time 
reversibility), the spectral statistics are conjectured to be of random-matrix type from the 
Gaufiian Orthogonal Ensemble [1, 2]. If the motion is generated by a Hamiltonian of a 
half-integer spin particle, then instead the statistics of the Gaufiian Symplectic Ensemble 
are expected [3, 4]. This is a strong form of universality, where a single exact expression 
for each spectral statistic is shared by all members of the universality class. 

Between the dynamical extremes of full chaos, and complete integrability, one can 
find quantum systems with a variety of intermediate spectral statistics. Systems with 
intermediate statistics exhibit properties such as linear level repulsion and an exponential 
decay in the probability for large spacings. However, as was pointed out in [5] they do not 
share precisely the same spectral statistics. 

Prominent examples that have been investigated in recent years include [6] the An- 
derson model at the metal-insulator transition point, [7, 8, 9] Aharonov-Bohm integrable 
billiards, [10, 11, 12] rectangular billiards, or other integrable dynamical system, perturbed 
by the addition of a point singularity (this example has often been called the Seha billiard) 
and [13, 14, 15] several other examples such as polygonal billiards with rational angles. 
(In [16] an alternative construction of the operator corresponding to the Seba billiard is 
given which exhibits Poissonian spectral statistics.) 

Some other situations in which intermediate spectral statistics appear include [17] the 
intermediate statistics of eigenphases of quantum maps, and [18] a one-dimensional model 



for intermediate statistics formed of a gas of energy levels interacting with a logarithmic 
potential. 

Spectral statistics of Seba billiards (see also [19, 20]) have been thoroughly investigated 
in [21, 22, 23, 24, 25, 26]. This model is particularly note- worthy since it constitutes an 
example of a member of a mini universality class of integrable systems perturbed by a 
rank-one singularity. Belonging to this class are quantum Neumann star graphs [27, 28], 
Seba billiards and quantum Neumann rose graphs (see section 2.2 below). 

As a measure of the spectral correlations we will mainly be using the pair correlation 
function R 2 {x) [29]. For a spectrum {A ra }, scaled so that the mean spacing is 1, we define 
R 2 (x)by 

J im m y^y^5(A„-A m ) =g(0)+ / g{x)R 2 {x)dx, (1) 

m=0n=0 J ~°° 

where g belongs to an appropriate class of test functions. The first term on the right-hand 
side of (1) comes from the diagonal terms m = n. By considering g to be approximately an 
indicator function we see that R 2 (x) is a measure of pairs of levels (regardless of ordering) 
that lie within a given distance of each other. Note that unlike the nearest-neighbour 
spacing statistic that is sometimes studied, R 2 (x) is not a probability density. 

The spectral statistics that are shared by Neumann star graphs and Seba billiards 
have been studied in some detail. The pair correlation function for small x was studied in 
[13, 22] and it was observed that it behaves as 

7T\/3 

R 2 (x) « — — x as x — > 0. (2) 

For the large x asymptotics, a full series expansion has been derived [27, 28]. Keeping 
terms up to x~ 12 the expansion reads 

„ , , 2 76 1088 9280 64000 n / 1 \ 

^) = 1 + ^2 + ^4-^6+^8- iio^ + O^n). asx^oo. (3) 

These behaviours are illustrated in figure 1, in which they are compared to a numerical 
calculation of the pair-correlation function. 

Our motivation is to analyse the spectral statistics of a system in the Seba class, with a 
symplectic symmetry. The most convenient model to use for this investigation are quantum 
graphs, since it is known that quantum star graphs with Neumann boundary conditions 
have statistics in the Seba class, and Bolte and Harrison [30, 31] successfully quantised 
quantum graphs with the Dirac operator, and showed that generically the statistics of the 
Gaufiian symplectic ensemble are found. 

The only obstacle to this programme is that the construction of [30] does not allow 
graphs with vertices of valency one, which is most of the vertices of a star graph. For this 
reason we re-attach the loose ends of the star graph to the central vertex, to form a graph 
which we call a rose graph (see figure 2). We will see (section 2.2) that despite this change 
in the graph topology important features of the spectral analysis for star graphs survive 
the transformation. 

The article is laid out in the following way. In section 2 we introduce the general scheme 
used to quantise the rose graph with the Dirac operator, in particular we derive a simple 
secular equation whose roots provide the spectrum of the graph with a similar structure 
to the well known secular equation of the Neumann star graph. Section 3 presents the 
derivation of the small and large parameter asymptotics of the pair correlation function 
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Figure 1: The pair correlation function numerically calculated for a B = 21 bond (green) 
and B = 101 bond (red) star graph with Neumann boundary conditions. Also plotted 
are the curves corresponding to the large and small values of the parameter, (3) and (2) 
respectively. The numerical plots are averaged over 100 realisations of the bond lengths, 
and 150 000 eigenvalues were computed. 



for the Dirac rose graph. Section 4 draws together the results and compares them with 
numerical calculations of the Dirac rose graph and a Neumann star graph. 

2 Quantisation of rose graphs with the Dirac operator 

A graph Q consists of a set of vertices V with pairs of vertices connected by bonds, as in 
figure 2. Two vertices u,v £ V are adjacent u ~ v if the pair (u, v) is in the set of bonds 
B, which may be associated with the set {1,2,... , B}. We will also use u ~ b and v ~ b to 
show that the bond b is connected to u and v. Each bond b is associated with an interval 
[0, Lb] where Lb is the length of b. Jzf = 2 J2b=i wu l denote twice the total length of the 
graph; a natural measure of the volume of the graph as each bond can be traversed in two 
directions. We will assume each Lb lies in an interval [1- (2B)- 1 , 1 + (2B)- 1 ] and that the 
set of bond lengths are rationally independent. On [0, Lb] we define a coordinate Xf, such 
that Xb = at the vertex o(b) and Xb = Lb at the vertex t(6); o(b) and t(b) are called the 
origin and terminus of b respectively. We will use x to denote a general coordinate on the 
graph when the particular bond is not significant. The number of bonds b with o(b) = v 
or t(6) = v is d v the degree of the vertex v. Note that a bond b with o(b) = t(b) = v, as 
seen in the rose graph, counts twice when determining the degree of v. It is sometimes 
convenient to have notation for the reverse of a bond; b is a bond with o(b) = t(b) and 
t(5) = o(6) (L- b = L b ). 

On an interval [0, Lb] the time independent Dirac equation reads, 

-i^ + m * = fi t. (4) 
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Figure 2: (a) A rose graph with 5 bonds; (b) a star graph with 5 bonds. 



where a and /3 satisfy the relations a 2 = j3 2 = I and a/3 + /3a = which define the 
Dirac algebra in one dimension. The dimension of a and j3 depends on the interpretation 
of the Dirac equation in one dimension. If one looks for the simplest faithful irreducible 
representation of the Dirac algebra a and (3 will be 2 x 2 matrices acting on two component 
spinors. On the other hand, if one regards the equation in one dimension as the restriction 
to a wire of a Dirac equation in three dimensions it is natural to expect four component 
spinors. It turns out that these two different approaches lead to the same scattering 
problem [30]. In order to impose time-reversal symmetry with two component spinors it 
is necessary to work with pairs of bonds connecting vertices, effectively reintroducing four 
components to the wave function for each adjacent pair of vertices. To simplify the current 
presentation we will only consider four component spinors from the outset and fix 
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Given the operator on the collection of intervals we must now specify a domain for 
which it is self-adjoint. To do this we fix appropriate matching conditions between the 
spinors on the set of intervals where they meet at the vertices of the Q. Such matching 
conditions were classified in [30] . In particular there is a natural generalisation of matching 
conditions that define a self-adjoint Laplace operator on the graph to conditions for a Dirac 
operator. 

For comparison the time-independent free Schrodinger equation on an interval reads, 



d 2 A 
dx 2 



(6) 



Matching conditions between the functions on the individual intervals where they meet at 
a vertex v are expressed as d v linear relations amongst the values of the functions ipb and 
their derivatives tfji evaluated at the end of each interval connected to v. For example, a 
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common choice of matching conditions for the Laplace operator are Neumann conditions; 
ijj is continuous at each vertex v; 



^j(O) = 7 for all b with 0(6) = v , (7) 
4> b (L b ) = 7 for all b with t(b) = v , (8) 



(Note that 7 is not a fixed constant but rather a placeholder for the value of the wave 
function at the vertex v whatever that happens to be) and the outgoing derivatives at v 
sum to zero, 

£ ^(0)- ti(L b ) = 0. (9) 

b\o(b)=v b\t(b)=v 

Fixing such matching conditions at each vertex of the graph ensures that the Laplace 
operator is self-adjoint. A full description of all matching conditions for which the operator 
is self-adjoint has been obtained in various forms [32, 33, 34], however for the current 
discussion an understanding of Neumann conditions is sufficient. 

For the Laplace operator it is clear that a wave function restricted to bond b can be 
written as a linear superposition of two plane waves, 

Mxb) = Q b e ikXb + R b e~ ikx " , (10) 

where k 2 = E. For simplicity we can consider the centre of the star v where v = o(b) 
for every b ~ v. Then Q = (Qi, ■ ■ ■ ,Qb) and R, defined similarly, are vectors of the 
coefficients of incoming and outgoing plane wave solutions at v. Using the matching 
conditions these can be related by a unitary matrix ctM (k) the vertex scattering matrix, 

R = a (v \k)Q. (11) 

With Neumann conditions at v the vertex scattering matrix takes a particularly simple 
form, see e.g. [35], 

>)_/ i [ib ^ c do) 

-{ £-1 if6 = c • ™ 

Components of cr^ are scattering amplitudes relating an incoming plane wave on bond c 
to an outgoing plane wave on bond b. Time-reversal symmetry implies that = (cr^) T . 

To quantise a graph with two component spinors one may treat the spinor matching 
conditions by analogy with those of the Laplace operator. The first component of the 
spinors evaluated at the ends of the bonds takes the place of the value of the scalar wave 
function of the Laplace operator while the derivative of the wave function evaluated at 
the ends of the intervals is replaced by the value of the second component of the spinor at 
the bond ends. With the four component spinors, necessary to incorporate time-reversal 
symmetry, the story is much the same. However, now a vector of the first and second 
components of the spinor is treated by analogy with values of the wave function and 
the third and fourth components replace values of the derivative. In the following we 
summarise the relevant results from [30] concerning the quantisation of a graph with the 
Dirac operator. 

If ij> = (-01, . . . iipk)^ is a 4-component spinor satisfying the Dirac equation (4) it can 
be expressed as a linear combination of four plane waves two traveling in each direction 
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along the bond. For E > m, l[> b is of the form 



4> 6 = Qb 
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(13) 



where k 2 = E 2 — m 2 , and Pb,Qb, Rb, Sb are constants of integration. It is convenient to 
define 



-iQb 

Pb 



and 



Tib 



iSb 

Rb 



(14) 



Then vectors of coefficients of the incoming and outgoing waves at the vertex of a star can 
be written as, 



Q = (Qi, ,Q 



B) 



and 



R — (R^, . . . , R^ 



(15) 



As in the scalar case, vectors of the incoming and outgoing coefficients can be related via 
a vertex scattering matrix; 

R = (J M (A;)Q . 



Time-reversal symmetry requires that 
where J is a block diagonal matrix 



(16) 
(17) 
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(18) 



If the vertex scattering matrix is also divided into 2x2 blocks (a^) bc , which relate pairs 
of incoming spinors on bond b to outgoing pairs on bond c, the symmetry condition (17) 
reads 

(o-M <) cb = det (V^) 6c ) (ia^^y 1 . (19) 
We can factor a block of the vertex scattering matrix so 

{a^) bc = X bc u hc , (20) 

where Ub c E SU(2), the symmetry condition is then X c b = Xb c and u c b = n^ 1 . In this 
notation spinor scattering on the graph is defined by a scalar scattering amplitudes Xb c 
which define the transition probability from bond b to c, just as in the case of scalar 
wave functions, and an additional rotation between the spinors when passing through 
the vertex, uj, c . In general these transition amplitudes and spinor rotations are obtained 
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from the matching conditions and will also satisfy the constraint that the whole scattering 
matrix er^ is unitary. One straightforward way to satisfy these symmetry relations also 
establishes a connection with the scattering matrices of the Schrodinger operator. A vertex 
scattering matrix of the form 



a 



Xj{v) ®I 2 } f tfW) 1 , (21) 



where X^ v ' is a symmetric unitary d v x d v matrix and = diag{m, . . . ,11^} with 
Ub £ SU(2) satisfies the symmetry condition (17). The d v matrices Ub can be thought of as 
defining a spinor rotation when leaving the vertex on bond b or the inverse rotation when 
entering the vertex from that bond. A symmetric unitary d v x d v matrix of transition 
amplitudes X^ 3 ' is precisely the form of vertex scattering matrix obtained for the Laplace 
operator with time-reversal symmetry (12). In fact, given matching conditions that make 
the Laplace operator self-adjoint with a d v x d v vertex scattering matrix X^ we can 
define analogous matching conditions which will make a Dirac operator self-adjoint for any 
given set of spinor rotations {u\, . . . , Ud v }- For instance matching conditions analogous to 
the Neumann conditions defined for the Laplace operator (which with a slight abuse of 
terminology we will also call Neumann conditions) are 

^o(6) v fc(0) = y for all b with 0(6) = v (22) 
u t(b)Vb{Lb) =y for all b with t(b) = v (23) 

where on each bond b we have such a solution \\> b = (ipbii ■ ■ ■ V'w)' 1 to (4) and we define 

= { ^ 2 (x) J ^ = { Mx) ) ■ ^ 

In (22)— (23) y is again a placeholder for a common value of the spinor at the vertex rather 
than a fixed constant vector. In addition the Neumann condition on the derivatives (9) 
becomes, 

^2 n o(6)W fe (0) - ^2 n t(6)W fe (L b ) = . (25) 

b\o(b)=v b\t(b)=v 

Together applying these conditions at all vertices of a graph defines a self-adjoint realiza- 
tion of the Dirac operator [30, section 5]. 

To each bond b of a graph we now see that there correspond a pair of 2 x 2 unitary 
matrices u Q (6) an d ^t(6) which rotate the spinor when entering the vertices at the origin 
and terminus of b respectively. At a vertex, during the transition from a bond bj to a 
bond bi the spinor is rotated according to the matrix u^b 6 SU(2) where 



Ubibi ■= «o(&0 (utfo) J , ( 26 ) 

It will be convenient to associate spin matrices to bonds rather than vertices, so for this 
reason we define 

Wb -=u o{ b) (27) 

and 

Wb ■= («t(6))~ u o(b) ■ (28) 
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We note that if b is the reversal of the bond b, 



w b = [ u t(b)) u o(b) 
= ( u o(b))~ u t(b) 

= w^\ (29) 

If the 2B matrices {u (&)> Ut(b)}beB are random and independently distributed, then the 
matrices {wb}beB are independently distributed. If the matrices {it (6) , Ut(b)}beB are inde- 
pendently distributed with Haar measure on SU(2), then the w b are Haar-distributed too, 
and the angles {9b}beB defined on [0, it] by 

Tr lib = Tr Wb =■ 2 cos 6b, (30) 

are identically independently distributed with a sine-squared distribution: 

2 f' x 

¥(9 b <x) = - sm 2 (f>d(j), O^x^vr. (31) 



o 
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The traces themselves have a semi-circle distribution: 



P(Tru; b <x) = — / \A 4 - t 2 dt, -2 ^ x 2. (32) 
2vr J_2 

Bolte and Harrison proved a trace formula for the eigenvalues of the Dirac operator on 
an arbitrary graph in terms of periodic orbits [30] . In order to state the trace formula, we 
need to define a few quantities. A periodic orbit p is a closed path on the graph, modulo 
cyclic shifts. So p can be written as a sequence of connected bonds p = (bi, 62, • • • , b n ) 
where t(bj) = o(bj + i) and t(b n ) = o(bo). We denote by the symbol V the family of 
all periodic orbits on the graph. It is possible that an orbit p € V can be written as 
a repetition of a shorter periodic orbit, in which case r p , the repetition number of p, is 
the maximal number of repetitions of a shorter orbit contained in p. A p will denote the 
product of all scattering amplitudes accumulated along the orbit. For the free Schrodinger 
operator this is simply be products of the elements of the vertex scattering matrices; 



While for the Dirac operator we replace cr^.^ with + ^ J . We will denote by d p the 



A _ /T (t(6i))_(t(6 2 )) (t(M) /OON 

A V - %6i a b 3 b 2 ■■■ a b 1 b n ■ \ 66 ) 

W with 

matrix product of all spin matrices accumulated on the orbit: 

dp ■= u bn b n ^ub n ^ lbn _ 2 ■ ■ ■ Ub 2bl u bl b n G SU(2). (34) 

Because of the cyclic invariance of the trace, we have 

Tid p = Tr(w bn w bn _ 1 ---w b2 w bl ). (35) 

£ p denotes the length of the orbit (i.e. the sum of lengths of bonds over which the orbit 
passes). 

Let h be a test function whose Fourier transform h(£) = h(k)e~ 2nlke dk is smooth 
and compactly supported. Then, in terms of quantities defined above, the trace formula 
reads: 

f; h(k n ) = *H(0) + 1^2- (^r) A P H(i p ), (36) 

n=l p^V P V 7 
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where 

H(£) := I h(k) cos k£dk. (37) 



In (36) recall that the quantity .if is twice the total length of the graph. 
2.1 Dirac rose graphs 

On the rose graph there are B intervals with both ends of each interval connected to 
the single central vertex, figure 2(b). Eigenspinors on the bonds have the form shown 
in equation (13). The Neumann matching conditions at the central vertex, which will 
determine the constants of integration, are, 

u o (b)Vb(0) = u t ( b )V 6 (£ 6 ) = y for all b, (38) 

which imposes continuity of the first pair of components of the spinor, and 

B B 

Y M o(fc)W fe (0) = ^2 u t{b)^b{L b )- (39) 

6=1 6=1 

Using the solution (13) and (14), condition (38) reads 



u o{b) (Q b + R 6 ) = u t[b) (Q b e ikLb + K b e~ ikLb ) = 



(40) 



+ m 
which gives 

«o(6) " u t(b) e ikL ^ Q b = - (u o(b) - u m e- ikLb ^j K b . (41) 



From (40), 



-JiL = Q 6 + R 6 , (42) 
and so, eliminating R;, using (41) we get 

Q b sinkL b = 2i / E + m u t(b)( u °Q>) ~ "t®^ ^oj)^ ( 43 ) 

Let us now, for simplicity, consider the case that L b k/ir Z for every bond b. Applying 
the boundary condition (39) gives 

B B 

-<(6)(Q6 - B*) = Y -<(b)(Q b e ikLb - R b e~ ikLb ), (44) 

6=1 6=1 

which simplifies to 

B 

J^2 (u o{b) - e ikLb u t{b) ) Q b = (45) 

6=1 



upon inserting (41). Substituting (43), we get 



1 B 1 



£ ^kT h ( u °® u m ~ eikLbh ) ( l2 " u ^( b f- ikLb ) y = o 



6=1 

B 



V . -. - {m + w^ 1 -2cosfcL fe I 2 )y = 0. (46) 



sin A;Lft 
=i 
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Since v)b + w b = Tiw^ = 2cos0&l2, the condition for k to be an eigenvalue becomes 

cos 6 b - cos kL b _ o ^ 
^ sin kL h 

Equation (47) is the secular equation for a Dirac rose graph. 

Let us now consider what happens if Lbk/ir = n E Z for some bond b. In that case, 
(43) may be re-written as 

= Q 6 sinA;X 6 = —j=== u -fo (u o(b) - (-l) n n t(fe) )n o( ) )) y 

= 2iVET^ U ^ ( '™ b " ( " 1)ni2)nt W U °(^- (48) 

Thus, if Wb 7^ ( — l) m l2, we find that y = 0, and consequently Q c = R c = for all other 
c ^ b, since the rational independence of bond lengths means that L c k/ir E" Z. From the 
boundary condition (39), we then would have Qb = R& = 0, and there are consequently 
no non-trivial solutions to (4). 

We therefore conclude that if all Wb ^ ±12, then the eigenvalues k n , n E Z are given 
by the condition (47), which is the central result of this paper. We also remark that if 
the matrices {^ (&) > u t(6) }&eB are chosen randomly independently with Haar measure, then 
almost-surely, Wb 7^ ±12- 

Kramer's degeneracy is evident in equation (46), since for any value of k satisfying (47), 
there is a 2-dimensional space of solutions for y . This leads to two linearly independent 
eigensolutions, so each eigenvalue occurs with multiplicity 2. 



2.2 Neumann rose graphs 

We can use (47) to derive the condition for k to be an eigenvalue of a rose graph quantised 
with the Laplace operator and Neumann boundary conditions. The key observation is 
that if we choose the SU(2) matrices to be 

u t(b) = ( n 1 ) f° r an bonds b, 



«o(6) - «t(b) _ I 1 

and choose m = 0, then the spectral problem given by (4) and (39) degenerates into a 
pair of independent copies of the Neumann eigenproblem for the Laplace operator. 

Since Wb = I2 for all b G B there will be eigenvalues of the form k = 2nir/Lb for n E Z. 
The rational independence of bond lengths implies that L c k/2ir E" Z for any other c 7^ b, 
so following the analysis of section 2.1 we get (c.f. equation (46)), 

B 1 

^ sinTr (2Is " 2coskL ^y = ( 49 ) 

c=l c 
c^b 

So y = 0, and thence Q c = R c = for all c ^ b. 

The remaining points of the spectrum come from the condition (47), which in this 
situation is expressible as 

y 1 - C ° 7 SkLb =yt a n(^-)=0. (50) 
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The equality here between the left-hand sides comes from the observation that 



1 " C0 ^ tan(£V iff^Z. (51) 



sin ip V 2 / 2ir 

What is particularly interesting is that condition (50) is exactly the eigenvalue condition 
for a Neumann star graph [28] with set of bond lengths {L^/2 : b £ B}. 

The spectrum {k n } of the Laplace operator on rose graphs with Neumann boundary 
conditions consists of points k = k n satisfying (50), together with points of the form 
k = 2mir/Lb with b £ B and m £ Z, which interlace the solutions to (50). Eigenfunctions 
corresponding to the latter class of eigenvalues are supported on a single bond b. 

3 Spectral statistics 

For a Dirac rose graph the eigenvalues are almost-surely the solutions k to the equation 
Z(fc) •= Z(k, L, 0) .= ^ gin(fcLfe) = 0. (52) 



Define 



so that 



costf-cosx 

z(x, 6>) := : , 53 

sin a; 



B 

Z{k) = Y J z{kL b ,9 b ). (54) 

6=1 

Let us note an alternative form for the function z(x, 9). From [36] (equations 1.421.3 
and 1.422.3) we have the pole expansions: 

^ oo 
k=l 



and 



Since 



we have 



1 ( 1 1 

coBecz = - + 2^ * . (56) 

fc=i 



1/1 &7T 1 kn , 

+ — + T-^oZ^ , (57) 



z 2 — k 2 ir 2 2 \z + kn 1 + k 2 ir 2 z — kir 1 + k 2 ir 
z(x, 8) = cos 9 cosec x — cot x 

= V ((-l)"cOBg-l)f— — ~ - m \ 2 ) 

rn — — rv~i ' 



(58) 



(Because of the regularisation, the series in (58) converges absolutely.) 

From this representation, we see that the structure of Z(k) on the real axis is a sequence 
of poles at the points k = mir/L), for m £ Z and b = 1, . . . , B. The only way that poles 
can become closely-spaced is if poles of z(kLb,9b) become close for different bonds 6. If 
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the set of bond-lengths {Li, . . . ,Lg} is linearly independent over Q, then the positions of 
the poles for different b become independent as k becomes large. 
By differentiating (52), we get 



B 

£ 

6=1 



Z\k) = YM l ~ C ° S ^ C :^ )>^ (59) 
sin [kLb) 



so that Z{k) is increasing between poles, and there is a unique zero of Z(k) between each 
consecutive pair of poles. 

These observations mean that spectral statistics of Dirac rose graphs fall into the class 
considered by Bogomolny et al. in [22]. By averaging over the random phases {9b}b£B, 
their analysis can be used — with only slight modifications — to derive an expression for the 
averaged pair-correlation function R2(x) (see equation (142), loc. cit.). In the following 
two sections we derive an approximation to B,2(x) for small x, by following a method 
developed in [22] and other places; and for large x, by a different method based on the 
trace formula (36). 

3.1 Small x behaviour of Rz(x) 

We shall use the method from [13, 22, 5] (see also [37]) that begins with the observation 
that statistics of small spacings of zeros of a random meromorphic function, real with 
poles on the real axis, are approximated by statistics of the zeros of a function with three 
randomly-distributed poles: 

A\ A 2 A 3 , , 

+ ~j + ~j = 0, (60) 



k — ei k — e2 k — e 3 

where Ai,...,A 3 are random amplitudes, and ei,...,es are random points. We can 
assume that the ei, . . . , e% are distributed over a wide interval of the real axis, since only 
close spacings will contribute to the approximation for statistics of small zeros spacings. 

From (52) and (58), the amplitudes Aj are given by the quantities (— l) mj cos8j — 1 
where mj G Z and 6j are independent identically distributed random angles according to 
(31). More explicitly, Aj is a linear transformation of cos8j, so it follows from (32) that 
the probability density of each Aj is 

2 



7T 



Vy(2-y)i[ ,2](y). (61) 



We follow the method of [22] which was used for a similar situation in which all Aj 
are equal to 1. We begin by re-arranging (60), to find that the solutions are given by the 
quadratic equation: 

(Ai+A2+A 3 )k 2 -(Ai(e2+e 3 )+A2(ei + e 3 ) + A 3 (ei + e2))k+Aie2e3+A 2 eie3 + A 3 eie2 = 0. 

(62) 

Denoting the solutions by ki t 2, we have 

, Ai(e 2 + e 3 ) + A 2 (ei + e 3 ) + A 3 {ei + e 2 ) ± VV 

h > 2 = 2(Ai + A 2 + A 3 ) ' (63) 

where 



V = {Ai{e2 + e 3 ) + A2{ei + e 3 ) + A 3 {ei + e2)Y-A{Aie2e 3 + A2eie 3 + A 3 eie2){Ai+A2 + A 



■'>) 



(64) 
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If L is any real number, then the translation (ei, e 2 , e^) 1— > (ei + L, e 2 + + £) shifts /ci 
and fc 2 by an amount L. It is therefore convenient to shift to a set of coordinates in which 

Ai{e 2 + e 3 ) + A 2 (ei + e 3 ) + A 3 (ei + e 2 ) = 0. (65) 

Then (64) becomes 

V = -4(A ie2 e 3 + A 2 eie 3 + A 3 eie 2 )(Ai + A 2 + A3) 

... 4 ^ A 2 (A 2 + A 3 ) 2 1 Ax (Ax + A 3 ) 2 2A!A 2 \ 

= 4(A X + A 2 + A 3 ) ^ ^ i + A2 ei + Ai + A2 e 2 + I7T ^e 1 e 2 j . (66) 

To get the leading contribution to the pair-correlation function, we average over the posi- 
tions ei and e 2 , and the random amplitudes Ai, . . . , A3: 



1 (- /'OO /'OO 

K(R 2 (x))^-eI / / <5(x - Ak) deide 2 

2 I J—00J—00 



(67) 



where 

2\/P 

Afc = k% — k\ 



2(Ax + A 2 + A 3 ) ' 

and T> is given by (66). To perform the integral in (67), we switch to a system of polar 
coordinates. Let 



/ Ai + A 2 / Ai + A 2 . , 

ei= V ^ 2 +^) pcosy? and e2= v ^i+^r iny - (69) 

Then we get 

,2tt roo Ai + ^ 



v / AiA 2 (A 1 +A 3 )(A 2 + A 3 ) 



o x 1 H cos <p> sm 0? dpdo9>, 71) 

V VAi + A 2 + A 3 V V / (A 1 +A 3 )(A 2 + A3) ' 1 1 rf 

which becomes, after a further substitution, 

ifiVB^- 1 !?/ Ai + A 2 A! + A 2 + A 3 

E(i? 2 (x)) « -E<^ 



2 I a/AiA 2 (Ai + A 3 )(A 2 + A 3 ) 4 

1+ /t— : — , = I =cos99siny dip p5(x-p)dp\. (71) 
'0 V V / (A 1 + A 3 )(A 2 + A 3 ) / 7 J 

We use the known integral formula 

'** dp 2vr 



o 1 + a cos p sin <p y 7 ! - a 2 /4 ' 



for |a| < 2, (72) 



to get 







27T / 2 V / ATA 2 " X" 1 / A X A 2 \~ 1/2 

1 -1 cos tp sin <p\ d<p = 2,tt I 1 

' ' ^l-t- ^3A^2 -t- ^3J/ 

(73) 



V / (ATTA3KA 2 TA37 r V r " V (A 1 +A 3 )(A 2 + A 

27r A /(Ai + A 3 )(A 2 +A 3 ) 



a' /2 VA! + A 2 + A 3 
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Substituting into (71) we arrive to 



As we are assuming that the Aj are identically distributed, we may symmetrize (74) and 
get 

, . , N ttx ( (Ai + A 2 + A 3 ) 3 / 2 ] . . 

E(fl2(I)) 53 y e { -limr }• (75> 

This coincides with equation (155) from [22], where it was stated but not explicitly derived. 

We have to evaluate (75) in the case that A\, . . . ,A 3 are independent and identically 
distributed with probability density (61). To do that, it is helpful to observe that since 
the Aj are identically distributed, we have 

E<J (Ai+A 2 + A 3 )^ \ = ( At 1 

VAxA 2 A 3 J I Al /2 Al /2 VAx + A 2 + A 3 



■,1/2,1/2 1 

+ 6E< jdi_^ — 1 y (76) 

1 ^Ax +A 2 + A 3 { 



Also, since Aj > almost-surely, we also have the identity 

1 1 



yMl + A 2 + A 3 ^ 
By Fubini's theorem, we have, from (76), 



/CO 
e - iAl+A2+A3)x 2 ^ (?7) 
-oo 



E 



(yLl + / 2 + ^ 3)3/2 ] = A r e|aP ^A-vv^+^n dx 



3 Z" 00 

-= / / 3/2 (x)/_ 1/2 (x) 2 + 2/ 1/2 (x) 2 /_ 1/2 (x)dx, 

V 71 " J-oo 



(78) 



where 



I„(a;) :=E{^e- Ai:c2 } = lj\»+V\2-y) 1 / 2 e-y x2 dy. (79) 
From [36] formula 3.383.1, we find that 

W= 7= pv ^;/ iFi(^ + |;^ + 3;-2x 2 ). (80) 

The asymptotic behaviour of the confluent hyper geometric function allows us to conclude 
that 

4/2(x)/-i/ 2 (x) 2 + 2/ 1/2 (x) 2 /_ 1/2 (x) = 0(4o ) asx^ioo, (81) 



so that the integral in (78) converges quickly, and may be accurately evaluated using 
numerical integration techniques. So doing, we arrive at the following result: 

7TC 

E(R 2 (x)) « — x as x -»• 0, (82) 
6 
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where 

c « 6.781 . . . (83) 

The numerical value for the slope in (82) is approximately 3.550 . . ., which is slightly larger 
than the corresponding slope for star graphs quantised with the Laplace operator (2). 

In order to verify the calculation of the constant c we have performed a Monte-Carlo 
simulation. By randomly generating 10 6 realisations of the random variable 

(A 1+ A 2 + A 3 )^ m 

we found a sample average of 6.785 with a standard deviation of mean of 3.296 x 10~ 3 . 
3.2 Large x behaviour of R2(x) 

We base our study of the pair-correlation function at large x on the trace formula. Starting 
from a smoothed version of (1) and applying a standard argument, using the trace formula, 
we get that 

. N N 

felArEE5( A n-A w ) 



iv-j-oo AN 

m=0 n=0 



■A-nAq^pEq ( Tr dp Tr dq \ _ / \ f £p \ A ( &p 



+ 9 



P,geV p i \ / \ \ / \ 

/•oo 

9(0) + 9(r)K 2 (T)dr, (85) 



where i^2( T ) is the distribution defined by 



1 v - ^ ApAqlplq i Tr dp Tr dq \ i . tp 



- i^^j s i w - ij (86> 

The quantity (86) is sometimes called the pair- correlation form factor and its behaviour 
as r 1 will determine the large x behaviour of R2(x). 

The leading-order behaviour of the form factor can be determined using Berry's di- 
agonal argument [38], and was investigated for generic graphs quantised with the Dirac 
operator in [30]. 

Let us denote by £ the set of all possible lengths of periodic orbits. Orbits with exactly 
the same length are said to belong to the same degeneracy class. We sort the sums in (86) 
by degeneracy class: 



1 X^&*(\ i A » T rdr- 



ee£ v 7 y P eV p 



2 



If we define 



2B ^ „o/ ^ ^pTrd p x 2 



r 



p 



where Vt is the set of periodic orbits of topological length t steps, then K(t, B) — > K2(t) 
weakly as B — > oo provided that t/2B — > r as B — > oo. (One proves this by integrating 
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(88) and (87) against a test function localised at r and noting that the bounds that we 
imposed on the bond lengths imply that the delta functions corresponding to orbits with 
different topological lengths are supported on disjoint intervals.) 

Let us order the sum over £ according to the number of distinct bonds to which the 
sum is confined, and average with respect to the random spin matrices. Then 

B 

E(K 2 (t))= ton ^KjfaB), (89) 

B— >co * — ' 

t/2B^r i=l 

where 

I restricted p€Vt ^ 

to j bonds e p =£ 

Our approximation to the averaged form factor will be based on selecting from the sums in 
(90) only those orbits which back-scatter the maximum number of times in their degener- 
acy class. As back-scattering is increasingly more strongly weighted as B — > oo for a rose 
graph, this may be expected to give a good approximation to the form factor. This type 
of approximation was considered for Neumann star graphs in [27], where it was compared 
to an exact expansion of the form factor, and was shown to reproduce exactly the first 
four terms of the Maclaurin expansion of the form factor. 

Let us first consider the special case j = 1 of periodic orbits confined to a single bond 
b of the graph. We will (unlike in [27]) need also to consider the parity of t. 

In the case that t is even, we can back-scatter t times by bouncing 1 back-and-forth 
along a single bond, so the contributing orbits are t/2 repetitions of bb (denoting by b the 
reversal of b) with B choices for the bond b. For such an orbit p the repetition number 
is r p = t/2, the stability amplitude is A p = (1 — 1/B) and d p = I2. Approximating the 
length of the orbit by t (as each bond length approaches 1 as B — > 00) we get 

lim Ki(t,B) « ton B—rt 2 f 1 - — ) -=■ ( - ) . (91) 
B^oo U ' 7 B-,00 Jz^ 2 V B t 2 \2 v ' 

t/2B^r t/2B^r x ' x ' 

As ££ -> 2B and t ~ 2Bt, we get 

/ 1 ^ iBr 
lim Ki(t,B) « lim 2 1 

B^oo B->oo V B 

t/2B^r v 

= 2e" 4r , (92) 

for t even. 

For t odd, we can back-scatter t — 1 times by bouncing back-and-forth and the final 
bond passed over is the same as the penultimate one; for each of B bonds, there are two 
orbits with maximal back-scattering, 

bbbb---bbb and bbbb---bbb. (93) 

For these two orbits, respectively Tr d p = Tr Wb and Trd p = Trw^ = Tr Wb, referring to 
(35) and (29). Both orbits have 

1 v- 1 1 



^=1^-1 5 («) 



1 We will refer to a transition of the form bb or bb as a bounce. 
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and £ p ~ t and r p = 1. Therefore, for t odd, we have 



t/IB^r t/2B 



1 / 1 \ 4 ^"^ 

= 2E(Trw 6 )V 2 e- 4T . (95) 
For a random matrix to £ SU(2) with Haar measure, E(Tru;) 2 = 1, so we get 

lim K x {t,B) « 2rV 4r . (96) 

B— >-oo 

Since odd and even values of £ have relative density 1/2, we find the total contribution for 
orbits confined to a single bond is the average of the two cases in (96) and (92): 

(l + r 2 )e- 4r . (97) 

We next will consider the case j = 2. This is more typical of the general case. We 
now consider orbits confined to two bonds, which we will denote by a and b. To maximise 
back-scattering, we bounce t% times on bond a and t2 times on bond b, so that t\ + £2 = t. 
In this way we can achieve t — 2 back-scatterings (since two transitions are necessarily 
between different bonds). 

For t\ and t% there are four possibilities:- 



t\ odd, t2 even (-g— cases) 
t\ even, ti odd (^i cases) 



t odd, 



2 

• ti even, £2 even (^f 2 - cases) 1 

jjj. jj )t \ r =^ * even. 

• t\ odd, E2 odd (2 cases) J 

Since odd and even values of t have relative density 1/2, each possibility is weighted 
approximately i/4 in its contribution to the form factor as t becomes large. 

In the first of the four possibilities listed above there are four members of the degeneracy 
class, which we may label as: 



aaaa ■ 


■ aaabb ■ 


■bb 


aaaa ■ 


■ aaabb ■ 


■bb 


aaaa ■ 


■ aaabb ■ 


■bb 


aaaa ■ 


■ aaabb ■ 


■bb 



The values of Tt d p corresponding to these four orbits reduce to (respectively): 

Trw a , 

Tru> s = Tr w a , , > 

(98) 

Trw a , 

Tr w a = Tiw a . 

Similarly, to the second of the four possibilities, the we find four members of the degeneracy 
class, each with Trd p = Trw D - 
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In the third possibilities, each bond is paired with its reversal an equal number of 
times, and the four members of the degeneracy class have Tr d p = TYI2 = 2. 
Finally, for the final possibility, the values of Tr d p are (respectively): 



Trw a w b , 
TvwaW b , 
Tvw a w- b , 
Trwawi. 



(99) 



This careful categorisation of each case reveals that the value of number of SU(2) matrices 
appearing in the expressions for TrcL depends only on the number of odd values of kj. 
In all cases, we have: 



1 v- 2 1 



A p=[b- 1 ) ( 10 °) 

i P ~ t, (101) 
r p = 1, (102) 

and there are ~ choices for the pair (a, 6), since the B(B — 1) free choices of 

bonds would lead to cyclic permutations, which should be counted only once in the trace 
formula. 

Putting the ingredients together, we have 
Jim K,{t,B)~ Jim ^^(l4) 2, " 2 i J |i( E (4^(^» f) +E (4^( T Vm^) 

t/2B-+r t/2B-rr 

+ 8 2 + E ((Tr w a w b + wg,w b + w a w b + w a wi) 2 ) "\ . (103) 

As before, E{(Tiw a ) 2 ) = E((Trw b ) 2 ) = 1, and we find that 

E ((Trw a w b + w a w b + w a w- b + WaW- b ) 2 ) = 4, (104) 
as a special case of equation (111) below. Therefore, 

8r 3/ lN 4rfl 



lim K 2 (t,B)& lim — 1-- (16 + 16 + 64 + 4) 

S->-oo B-i-oo 64 V 

t/2B^r v y 

= f r 3 e^. (105) 

We now turn to the general case of orbits confined to j bonds. An orbit with maximal 
back-scattering will bounce along bond d£ a total of ti times, I = 1, . . . , j, in such a way 
that 

ti H h tj = t. (106) 

In this way, t — j back-scatterings are achieved. 

The weighting that each possibility receives in its contribution to the form factor has 
two components. The first component is based on the relative density of occurrences of 
odd or even tis. Since there are are 2 J possible choices for ti to be odd or even, this factor 
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is simply 1 j2 3 . A second factor comes from the number of ways to decompose t in the 
form (106). The number of such decompositions is the number of ways to choose j ' — 1 
numbers (the transition points: ti, t\ + *2, et cetera) from a total of t — 1 possible ones. 
Although these two factors are not independent, since the odd and even choices for tg are 
evenly distributed, the weight of each possibility is approximately their product: 

1 (t-\\ P' 1 , s 

(107) 



z>(j-iy. 

Within each possibility it is the number of odd tjS that determines the value of the 
trace factor; after cancellations we are left with 

Tr d p = Tr (w% ■■■w%? r ), (108) 



where there are precisely r indices . . . ,i r } for which ti t is odd, and each ati e = ±1. 
There are 2 3 members of each degeneracy class, and 2 r ways that = ±1. So we will 
need to calculate 



E 



{(^ £ Tr J. (109) 

Since every combination of = ±1 occurs in the sum, and w ai + = Tr(w ai we 
have 

w:il--- W :r r =(Tiw ai )---(Trw ar )h. (110) 

a i( =±l 

Thus, 

E {f 2J " r E Tr(^ •••<:)) } =V +1 - r E{(Trw ai ) 2 -.-(Tvw ar f} 



a if =±l 



A j+l ~ r . (Ill) 



Finally, we note that there are r ) combinations of ways that there can be r odd indices 
out of a total of j. 



In all cases, we have: 



I \'-> i 



l p ~ t and (113) 
r p = 1. (114) 

1 B\ Bi 



The number of choices for the j bonds on which the orbits are confined is 

j(B-j)\ j 
(the factor 1/j is to account for cyclic invariance). We get 

t/2B^r V 7 J lj=0 V 7 J 

= t^'+V 4 ^. (115) 
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We remark that upon substituting j = 2 into (115) we recover (105), as expected. 
Substituting (97) and (115) into (89) we get 



5^' 



E(K 2 (t)) « (1 + r 2 )e" 4r + ^r- 7 ' +1 e- 4r 



i=2 J 



= (1 + r 2 )e" 4T + re" 4T (e 5r - 1 - 5r) 

= (l-r-4r 2 )e- 4r + re T . (116) 

Expanding (116) as a Maclaurin series, we get 

1 3 

E(# 2 (r)) « 1 - 4r + 9r 2 - -V + 0(r 4 ), (117) 

6 

for small values of r. 

It follows from the fact that the form-factor and the pair-correlation function are 
related via the Fourier transform that the small r asymptotics (116) determine the large x 
behaviour of the averaged pair-correlation function. Namely, we have [39, page 102] that 
if k(r) is even and 



and 



then 



fc(r)~l + ^a,r fc (118) 

k=l 

/•oo 

l-fc(r)= / (1 - r(x))e 2nixT dx, (119) 



-w~^* e (f;(gr^ f . 



(120) 

V '— I r" ■ • 1 

k=l 

Applying this to the form-factor approximation (117), we get the following approximation 
to the pair-correlation function: 

2 13 / 1 \ 

w'^-^+oy, (mi 

for large values of x. 

4 Conclusions 

We have investigated the spectral statistics of a quantum graph quantised with the Dirac 
operator, for which the Schrodinger-operator-quantised counterpart (Laplace operator) 
has intermediate spectral statistics. The shape of the graph has a single central vertex, 
and all bonds are connected at both ends to that vertex. We call this graph a rose. 

We have shown that the generic condition for k to be an eigenvalue is that it satisfies 
the non-linear equation 

^cosg b -cosfcL b 

^ sin kL b 

6=1 

where the angles Q\, are determined by a set of B matrices from SU(2) that rotate the 
spinor during its passage along the bonds. Each eigenvalue occurs with multiplicity 2 
(Kramer's degeneracy). 
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We have investigated the behaviour of the spectral pair correlation function R.2(x) 
in the limit B — > oo for large and small values of x, when it is averaged over random 
realisations of the SU(2) matrices chosen with Haar measure. 

For small values of x we have found 

7TC 

E(R 2 (x)) » —x asx^O, (123) 

where c ps 6.781 . . . 

For large values of x we have determined 

2 13 / 1 \ 

E(fl2W) . 1+ ___ +0 y . d24) 

The behaviour of the pair correlation function is different to the pair correlation func- 
tion for star graphs quantised with the Laplace operator (see also figure 4 below for a 
numerical comparison). It is too early to conjecture that the behaviour that we found 
is universal for systems with intermediate statistics and a symplectic symmetry, but this 
possibility merits further investigation. 

It is interesting to note that while we considered elements Wb chosen with Haar measure 
in SU(2) the large parameter asymptotic of R2(x) would be the same if Wb is chosen from 
any irreducible representation of a subgroup r C SU(2). The reason for this was seen in 
(111) where to evaluate the asymptotic it was only necessary to know that E{(Trw) 2 } = 1. 
As Tr w is the character of an element of the subgroup when the average is carried out 
over r for an irreducible representation of r the result must still be one by the character 
orthogonality relations. So, for example, choosing the spin transformations Wb on the rose 
from the finite subgroup 

r = {±l,±ia x ,±ia y ,±ia z } , (125) 

where Oj is a Pauli matrix, will not change the large x asymptotic of i?2(x). The small 
parameter asymptotic, in contrast, depends on the distribution of the Aj given in (61) 
which will vary if the spin transformations are chosen from an irreducible representation 
of a subgroup. 



4.1 Some numerics 

In order to numerically verify the small and large x behaviour of the pair correlation func- 
tion found in the previous section, we performed numerical calculations of the eigenvalues 
of rose graphs quantised with the Dirac operator, and calculated the empirical pair corre- 
lation statistic. To implement the averaging over the random choice of SU(2) matrices at 
the graph vertex, we performed the calculations 100 times with random Haar-distributed 
matrices and averaged the results. In each realisation, 150 000 eigenvalues were calculated. 

In figure 3 we compare the numerically-calculated pair correlation function with the 
predictions of (123) and (124), and find a good agreement in the range of validity. 

The agreement gets better as the number of bonds increases, as is to be expected since 
our analytical calculations relate to the limit B — > oo. This improvement is demonstrated 
in figure 6 in which a comparison is made between the pair-correlation function for 21, 61 
and 101 bond graphs, and the large x prediction of (124). A clear increase in adherence 
to the prediction is displayed as the number of bonds increases. 

One may wonder if averaging over both the random SU(2) matrices and bond lengths 
will lead to a different pair-correlation function. Figure 5 compares such a calculated 
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empirical curve with that for corresponding graph with averaging over the SU(2) matrices 
only, for a B = 101 bond Dirac rose graph. There is no noticeable difference in the curves 
so obtained. 

It is of interest to compare the pair correlation function for the Neumann star graphs 
with the Laplace operator, and the Dirac rose graphs. This is done in figure 4 for graphs 
with 101 bonds in both cases. Qualitatively the curves are similar in appearance, but 
there is some noticeable difference, which can already be explained analytically around 
the point x = (compare (123) and (2)). 
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Figure 3: The pair-correlation function numerically calculated for a B = 21 bond (green) 
and a B = 101 bond (red) Dirac rose graph. Also plotted are the curves corresponding to 
the large and small values of the parameter, (124) and (123) respectively. The numerical 
plots are averaged over 100 realisations of the SU(2) matrices, and 150 000 eigenvalues 
were computed. 
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